1 Importations et transformation des données

1.1 Importations

1.2 Sélection des données

## [1] <NA> 1    2    3    4   
## Levels: 1 2 3 4

On récupère les données suivantes :

  • obsdata, segdateet predata
  • distdata, une jointure entre predata et obsdata sur le segment

1.3 Centrage-réduction

On centre et réduit les covariables (présentes dans distdata, segdata et predata).

##      CHL_4w_mea CHL_4w_sd SST_4w_mea SST_4w_sd POC_4w_mea    depth distCoast
## mean  1.3016453 0.4070661  15.865554 0.8412628  211.17688 51.14547  38510.47
## sd    0.8536337 0.3629038   3.251786 0.4386712   94.86176 36.32341  30059.10
##       dist200   slopeP
## mean 88112.13 28913.18
## sd   32477.64 27947.08

1.4 Longitude, latitude en X et Y en lambert93

On recalcule X et Y en lambert 93 à cause d’un souci dans les données.

predata

distdata

gridata

1.5 Dates des sessions

## # A tibble: 4 x 2
## # Groups:   session [4]
##   session dates_session                                                         
##   <fct>   <chr>                                                                 
## 1 1       2019-02-12 + 2019-02-13 + 2019-02-14 + 2019-02-15 + 2019-02-26 + 2019~
## 2 2       2019-05-30 + 2019-05-31 + 2019-06-01 + 2019-06-02 + 2019-06-03        
## 3 3       2019-07-31 + 2019-08-01 + 2019-08-02 + 2019-08-05 + 2019-08-07 + 2019~
## 4 4       2019-10-25 + 2019-10-26 + 2019-10-27 + 2019-10-28 + 2019-11-17 + 2019~

2 Fonction de détection

2.1 Sélection

On cherche à savoir quelle est la meilleure fonction de détection.

## [1] "Itération 1 sur 8 terminée."
## [1] "Itération 2 sur 8 terminée."
## [1] "Itération 3 sur 8 terminée."
## [1] "Itération 4 sur 8 terminée."
## Error : 
## gosolnp-->Could not find a feasible starting point...exiting
## 
## [1] "Itération 5 sur 8 terminée."
## [1] "Itération 6 sur 8 terminée."
## [1] "Itération 7 sur 8 terminée."
## [1] "Itération 8 sur 8 terminée."
##   observerId seaState              formula key       AIC
## 1          0        0                   ~1  hn -293.4346
## 2          0        1            ~seaState  hn -290.1428
## 3          1        0          ~observerId  hn -288.9999
## 4          1        1 ~observerId+seaState  hn -288.8973
## 5          0        0                   ~1  hr -297.3398
## 6          0        1            ~seaState  hr -297.3362
## 7          1        0          ~observerId  hr -295.5575
## 8          1        1 ~observerId+seaState  hr -295.8082

On choisit la fonction de détection avec l’AIC le plus faible. Ici, il s’agit de la fonction avec la formule ~1 et la key hr. A cause d’une erreur pour ajuster cette fonction seulement, on préferera prendre la fonction avec la formule ~seaState et la key hr, qui a un AIC très très proche. On sélectionne donc la meilleure fonction de détection : detfc.sea.hr

## 
## Summary for distance analysis 
## Number of observations :  96 
## Distance range         :  0  -  0.282 
## 
## Model : Hazard-rate key function 
## AIC   : -297.3362 
## 
## Detection function parameters
## Scale coefficient(s):  
##               estimate        se
## (Intercept) -2.2314812 0.1847034
## seaState     0.1731484 0.1239956
## 
## Shape coefficient(s):  
##             estimate       se
## (Intercept) 1.473509 0.242613
## 
##                        Estimate          SE         CV
## Average p             0.5468238  0.04173954 0.07633087
## N in covered region 175.5593056 18.07215909 0.10294048

2.2 Résultats numériques et graphiques

2.2.1 Avec plot.ds

2.2.2 Manuellement, avec ggplot

On calcule manuellement les courbes de détection en fonction de l’état de la mer.

# On récupère l'ensemble des paramètres du modèle d'abord
detfc_par <- detfc.sea.hr$ddf$par

# scale.value : fonction de mise à l'échelle des paramètres du modèle (matrice et exponentielle)
scale.value <- function (param, z){
    exp(as.matrix(z) %*% param) }

# key.fct.hz : On entre une distance, un sigma (scale), un b (shape) et on obtient une probabilité de détection pour une hazard-rate
key.fct.hr <- function (distance, sigma, b){
    return(1 - exp(-(distance / sigma) ^ (-b))) }

# On créé un vecteur des distances de 0 à 0.3 km
dist <- seq(from = 0, to = 0.3, length.out = 96)

# On récupère le paramètre shape (b) et on le met à l'échelle : 100 fois exp(b)
b <- scale.value(detfc_par["V1"], matrix(1, nrow = 96, 1))

# On va ensuite récupérer le paramètre scale (sigma)
# C'est plus compliqué ici, car sigma varie selon les covariables de détection, donc seaState ici.

# sigma.0 : intercept
sigma0 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState = 0
                                            z = matrix(c(rep(1, 96), rep(0, 96)), ncol = 2))

# sigma.1 : intercept + seaState toujours à 1
sigma1 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState = 1
                                            z = matrix(rep(1, 96*2), ncol = 2))

# sigma.2 : intercept + seaState toujours à 2
sigma2 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState constant à 2
                                            z = matrix(c(rep(1, 96), rep(2, 96)), ncol = 2))

# sigma.3 : intercept + seaState toujours à 3
sigma3 <- scale.value(detfc_par[2:3], # 2 paramètres : sigma intercept + sigma seaState constant à 2
                                            z = matrix(c(rep(1, 96), rep(3, 96)), ncol = 2))


# A partir de tous ces paramètres et du vecteur des distances, on calcule la probabilité de détection
proba_detec0 <- key.fct.hr(distance = dist, sigma = sigma0, b = b)
proba_detec1 <- key.fct.hr(distance = dist, sigma = sigma1, b = b)
proba_detec2 <- key.fct.hr(distance = dist, sigma = sigma2, b = b)
proba_detec3 <- key.fct.hr(distance = dist, sigma = sigma3, b = b)

On va ensuite générer le graphique correspondant.

On enregistre ce graphique dans le dossier img en 2 versions, en anglais et en français, pour le poster.

## png 
##   2
## png 
##   2

Nous avons aussi tenté en mettant les modalités beaufort 0-1 ensemble et 2-3 ensemble, pour voir si nous avions toujours une meilleure détection avec une mer plus agitée.

3 Fonction de densité (covariables communes par session)

3.1 Analyse préliminaire : corrélation des covariables

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 25992 individuals, described by 9 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

3.2 Sélection des covariables

On utilise la fonction selec_dsm_aic_fwd qui va permettre de sélectionner de manière forward les covariables à inclure dans le modèle dsm.

Cette fonction prend en arguments segdata, obsdata, la fonction de détection et un vecteur de toutes les covariables à tester. C’est une fonction récursive, elle peut donc aussi prendre en argument un vecteur des covariables déjà sélectionnées.

On sélectionne les covariables sur l’ensemble des données, sans distinction de session. Ici, les AIC sont données pour availability = 1. L’algorithme pour les autres valeurs de disponibilité sélectionne les même covariables, seul l’AIC change.

## [1] "----- Sélection de la covariable numéro 1 -----"
## [1] "----- Sélection de la covariable numéro 2 -----"
## [1] "----- Sélection de la covariable numéro 3 -----"
## [1] "----- Sélection de la covariable numéro 4 -----"
## [1] "L'AIC ne diminue plus en ajoutant une 4ème covariable"
## $dsm.selec
## 
## Density surface model
## Response :  count 
## 
## Detection function :  Hazard-rate key function 
## 
## Formula:  count ~ s(SST_4w_mea) + s(X, Y) + s(CHL_4w_mea) + offset(off.set) 
## 
## Fitting engine: gam
## 
## 
## $aic
## [1] 1757.999
## 
## $formule.dsm
## [1] "count ~ s(SST_4w_mea) + s(X, Y) + s(POC_4w_mea)"
## 
## $suivi
##   iter                                         formule      aic
## 1    1                           count ~ s(SST_4w_mea) 1842.067
## 2    2                 count ~ s(SST_4w_mea) + s(X, Y) 1763.128
## 3    3 count ~ s(SST_4w_mea) + s(X, Y) + s(CHL_4w_mea) 1757.999

3.3 Ajustement de la fonction de densité

dsm pour \(availability\) dépendante de on-shelf et off-shelf : On note “on-shelf” quand la profondeur est inférieure à 150m, et “off-shelf” si la profondeur est supérieure à 150m.

\[availability_{off-shelf}=0,1357617\] \[availability_{on-shelf}=0,6332016\]

4 Prédiction de l’abondance et de la répartition

On créé deux tableaux de données predata_tmp dont les valeurs correspondent aux valeurs de predata pour chaque session.

4.1 Abondances estimées selon la disponibilité et la session

##    Session      Disponibilité       Min Estimation       Max
## 1:       2                  1 2948.3582   4441.895  5935.433
## 2:       3                  1  348.7289   1536.801  2724.874
## 3:       2               0.41 9170.5957  10664.133 12157.670
## 4:       3               0.41 2515.4844   3703.557  4891.629
## 5:       2 on-shelf/off-shelf 5452.3439   6945.881  8439.418
## 6:       3 on-shelf/off-shelf 1194.0925   2382.165  3570.237
##    Session      Disponibilité Estimation Plus.ou.moins
## 1:       2                  1   4441.895      1493.537
## 2:       3                  1   1536.801      1188.072
## 3:       2               0.41  10664.133      4266.205
## 4:       3               0.41   3703.557      3566.324
## 5:       2 on-shelf/off-shelf   6945.881      2531.880
## 6:       3 on-shelf/off-shelf   2382.165      2074.952

4.2 Cartes statiques (poster)

On créé maintenant toutes les cartes statiques. Ce sont celles que l’on retrouve sur le poster. Elles vont également être enregistrées dans le dossier img.

4.2.1 Session 2 et availability = 1

## png 
##   2

4.2.2 Session 3 et availability = 1

## png 
##   2

4.2.3 Session 2 et availability = 0.41

## png 
##   2

4.2.4 Session 3 et availability = 0.41

## png 
##   2

4.2.5 Session 2 et availability = off/on shelf

## png 
##   2

4.2.6 Session 3 et availability = off/on shelf

## png 
##   2

4.3 Cartes interactives (application shiny)

On obtient aussi une carte interactive par session, qui sont celles que l’on retrouve dans l’application.

4.3.1 Session 2

4.3.2 Session 3

5 Effets des covariables sélectionnées

5.1 Visualisation des splines correspondantes au dsm

On enregistre plusieurs images dans le dossier img : par session et biais de disponibilité, on a une image avec toutes les splines et une image avec X et Y en plus grand.

Voici par exemple l’image pour la sessio 2 (printemps), un biais de 1 et toutes les covariables.

5.2 Cartes des covariables SST et Chlorophylle

On récupère les predata avec les covariables non centrées-réduites.

5.3 Chlorophylle

5.4 SST

On retrouve toutes les cartes ci-dessous :

6 Exportation des résultats en Rdata

On exporte certains résultats en Rdata qui seront utilisés pour l’application shiny.